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Exploiting alternative representations for the Wigner lunction's dynamical equation, we develop 
a general strategy for effective numerical propagation. As an example, the split-operator approach 
is implemented and illustrated for a variety of systems. 
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The Wigner function was introduced at the dawn of quantum mechanics [1] resulting in a powerful tool to analyze 
the quantum-classical transition [2H1]- The Wigner function's dynamic equation, known as the Moyal's equation 
of motion OH], can be written as either a partial differential equation containing infinite-order derivatives, or an 
integro-differential equation. Hence, development of efficient numerical methods has been hindered by the equation's 
inherent complexity. 

Many different numerical methods were proposed based on solving the integral representation of the Moyal's equation 
[TlfTO] , reducing to an effective Boltzman transport equation ^TTJ [T2] , propagating Gaussian and coherent states [13l - 
[T5] . combining spectral methods with optimization techniques [TB], etc. 

Moyal's equation of motion is [T7] 
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where {{ , }} denotes the Moyal bracket. An explicit expansion of the bracket leads to 
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which can be expressed in terms of the generator of motion G as 



^^dw^^awix,p), (3) 

G^Hi^x-p,p+^X^-H(^x+ p,p - ^a) , (4) 
where the operators x,p,9,X satisfy the following commutator relations |18| 

[x,p] = o, [x,x\ = i, [p,e] = i, [xJ]^o. (5) 

Note that Eq. pi contains no direct reference to specific differential operators. 
The generator of motion for the Hamiltonian H = |^ -f U{x) reads 

G ^ —pX + U{x~h§/2)-U{x + he/2). (6) 

TO 

Moyal's equation ([T]) is recovered in the phase space representation x ~ p for which 

X ^ X, p^p, A = ~idx, = -idp. (7) 
Utilizing an alternative representation in the x — 9 space 

X = X, p = idg, X — —idx, 9 — 9, (8) 
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FIG. 1: The evolution of the Wiener function at times t=0 a.u. and t = 7.5 a.u. for the Morse potential U{x) = 20(1 - e"" "^). 
The initial state corresponds to the first exited state of the Harmonic oscillator. 



the equation of motion turns into 
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where the underlying Wigner function is recovered through a Fourier transform 

Wix,p) ^ J dee'P''B{x,e). (10) 

Moreover, there are two additional representations in terms of the A — 6* and the \ — p variables. The corresponding 
representations of the quantum state Z{X,6) and A{\,p) are obtained as 



Z{\,e) = J dxe-'^''B{x,e), 
A{X,p)^ [ dxdee''-P''-^'='>B{x,e). 
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The latter transformation is specially important and is hereafter denoted by J-', such that 



A{X,p) ^TB{x,e). (13) 

In contrast to Moyal's equation ([I]) in the phase space, Eq. ^ can be numerically propagated by a broad variety 
of numerical methods developed for the Schrodinger equation. In particular, the first- and second-order split-operator 
methods given by 
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can be effectively implemented by utilizing the Fast Fourier Transform from the x — 9 to X — p representation. For 
example, the first-order propagation takes the form 



B(t + At,x,e) 



(16) 



More details and a Python source code is provided in the Appendix. This method was employed to generate animations 
of the Wigner function propagation in the case of quartic [12], Morse |^, and Gaussian [5T] potentials. In particular. 
Fig. [T] shows the initial and final state of the Morse potential propagation. 

Conclusions. Exploiting the ability to express the time-evolution equation in alternative representations (where 
the phase space is a special case), we presented an effective and straightforward strategy for the Wigner function 
propagation. In fact, numerical evaluation is reduced to solving a Schrodinger-like partial differential equation. One of 
such methods is the spectral split-operator, employed to calculate the evolution for the Morse, quartic, and Gaussian 
potentials. This method can be generalized to non-unitary evolution. Moreover, a similar numerical approach is 
applicable to the spinorial relativistic Wigner function |22j. 
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APPENDIX 



The following is the Python implementation of the first-order split-operator. 

# Loading packages 

import numpy as np 

import scipy . fftpack as fftpack 

# Defining the potential 

def Potential (x): 

return 0.1*x**4 



# Specifying parameters in atomic units 



dis c r c t i z at io nX = 512 # Number of points in x 

d i scr e t i z at ionP = 512 # Number of points in p 

dt = 0.02 # time increment 

timcStcpsN = 180 # number of propagation steps 

mass =1. # particle 's mass 

amplitudcX = 9.0 # x range = [— amplitudeX , amplitudeX ] 

amplitudeP = 25.0 # p range = [— amplitudeP , amplitudeP ] 

hbar =1. # Planck constant 

# Defining x_vector , p_vector, theta_vector and lambda-vector 

# containing the respective range of values 

x_vector = \ 

np . linspace (—np . abs ( amplitudeX ) , np . abs ( amplitudeX *( 1 . — 2./ discretizationX ) ) , discretizationX ) 
p_vector = \ 

np . linspace (—np . abs ( amplitudeP ) , np . abs ( amplitudeP *( 1 . — 2./ discretizationP ) ) , discretizationP ) 

theta_vector = f ft pack . f f t s h i f t ( \ 

2 . * np .pi * f f t pack.fftfrcq( p_vector.size, p_vector[l]— p_vector[0] ) ) 

lambda_vector = fftpack . f ft shift ( \ 

2.*np.pi* fftpack . fft freq ( x_vector . size , x_vector [1] — x_vector [0] ) ) 

# Defining X,P, Lambda and Theta grids 

Tlicta , X = np.meshgrid ( theta_vector ,x_vcctor) 
P, Lambda = np.meshgrid (p_vector , lambda_vector ) 

# Defining the propagator factors 

potentialPropagatorFactor = fft pack . i f f t s h i f t ( \ 

np.exp( -lj*dt*(Potential(X-hbar*Theta/2.) - Potent ial (X+hbar*Theta / 2 .) ) ),axes = (l,) ) 

kineticPropagatorFactor = fft pack . i f f t s h i f t ( \ 

np . exp ( — Ij * Lambda*P*dt/mass ) , axes = (0,)) 

# Secifying the initial state 

W-init = np.exp( -(X-2)**2 - (P )**2) + Oj 

# Propagation 

W= np . array ( W_init ) 

for t in range ( 1 , timcStcpsN ) : 

# Transforming to the X— Theta representation 
B = fftpack . fft (W, overwrite_x=True , axis = l) 

# Applying the factor associated with the potential 
B *= potentialPropagatorFactor 



# Transforming to the Lambda— P representation 
Z = fftpack . fft (B, overwrite_x=True , axis=0) 
A= fft pack . i f f t (Z , overwrite_x=True , axis=l ) 
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# Applying the factor associated with the kinetic energy 
A*= kineticPropagatorFactor 

# returning to the X—P representation 

W = fftpack . iff t (A, overwrite_x=True , axis=0) 

W_final = W 

# Ploting the final state in the Python SAGE distribution (www.sagemath.org) 

matrix_plot( np.array(W_final.real) , cmap= ' hsv ' , color b ar=True , origin='lower ') 
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